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A detailed description and validation of a recently developed integration scheme is here 
reported for one- and two-dimensional reaction-diffusion models. As paradigmatic ex- 
amples of this class of partial differential equations the complex Ginzburg-Landau and 
the Fitzhugh-Nagumo equations have been analyzed. The novel algorithm has precision 
and stability comparable to those of pseudo-spectral codes, but it is more convenient to 
employ for systems with quite large linear extention L. As for finite-difference methods, 
the implementation of the present scheme requires only information about the local envi- 
roment and this allows to treat also system with very complicated boundary conditions. 

Keywords: Partial Differential Equations; Reaction-Diffusion Models; Integration 
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1. Introduction 

The study and development of efficient and accurate integration algorithms 
for reaction-diffusion partial differential equations (RDPDE) is a compelling 
task. Because reaction-diffusion equations are used to represent quite a large 
class of systems ranging from chemical, to biological and physical ones El. Within 
the field of physics deterministic reaction-diffusion models have been used to 
mimick front-propagation, surface- and interface-growths, turbulent behaviour. 



pattern formation, excitable media, etc. 

A paradigmatic example of RDPDE is the well-known Fitzhugh-Nagumo 
equation (FHNE) used to reproduce the onset of solitary waves, periodic wave 
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trains, circular and spiral waves in excitable media 

ut ~ DV^u + u{u — a){\ — u) ~ V ; vt = e{(3u — v) (1) 

where u — u{r, t) and v — v{t, t) are real fields and represent the activator and 
the inhibitor, respectively. The appearance of different dynamical behaviours 
and of the associated patterns is ruled by the values assumed by the real pa- 
rameters D, a, e and (3 and by the dimensionality of the system (see Refs. cli 
for a review). 

Another widely studied reaction diffusion model is represented by the com- 
plex Ginzburg-Landau equation (CGLE), that is particularly relevant for the 
description of the dynamics of spatially extended systems both from an ex- 
perimental and from a theoretical point of view. This because any spatially 
extended systems that undergoes a Hopf bifurcation from a stationary to an 
oscillatory state is described by the CGLE sufficiently close to the bifurcation 
point 0. The CGLE can be written as 

At = (1 + JCi)V^A + A + (1 - jc3)| (2) 



where ci and C3 are real positive numbers, while A{r,t) — p(r, exp [j'0(r, t)] 
is a complex field of amplitude p and phase V'- Here j is the imaginary unit 
with p = —1. Equation (|^) exhibits several different stationary and turbulent 
regimes and a quite rich literature has been devoted to the description of its 
phases B0. In the limit (01,03) (00,00) the CGLE reduces to the nonlinear 
Schrodinger equation, while the real Ginzburg-Landau equation is recovered in 
the opposite limit (01,03) (0,0) 0. 

In the present paper we describe in detail the implementation and the per- 
formances of a new integration scheme for RDPDE. In particular, we deal with 
the FHNE and the CGLE in one and two dimensions for different types of 
boundary conditions. In Section II the new integration algorithm is explained 
and its apphcation to the FHNE and to the CGLE is described in detail. The 
treatment of different boundary conditions (namely, periodic, no-fiux and fixed 
ones) is explained in Section III. In Sect. IV the precision and the performances 
of the present scheme are compared with those of usual pseudo-spectral codes 
for one- and two-dimensional systems. A brief final discussion is reported in 
Sect. V. 



2. The integration scheme 

The algorithm here reported is a modification of the well known "leap-frog" 
method applied to a spatially extended system, whose evolution is represented 
by a partial differential equation XPDE). The leap-frog algorithm was originally 
devised for Hamiltonian_systems 113 and heavily employed for classical molecular 
dynamics simulations lill. But it can be also implemented for any ordinary 
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differential equation, when the time evolution of a given dynamical variable 
z = z{t) can be written as 

z(i) = A z(<) + B z{t) (3) 

where A and B are two operators, generically non linear, that typically do not 
commute. Suppose now that one is able to integrate separately the two contri- 
butions associated to A and B in Eq. @, but not Eq. as a whole. Therefore 
one would like to rewrite the formal solution of @ z(t) = exp [t{A + B)]z{0) 
as a product of terms exp [afc^A] and exp [bktB]. A vast literature has been de- 
voted to the choice of the coefficients ak and bk that minimizes the errors done 
in rewriting the formal solutions in such an approximate way E3'E3. The best 
known approximation is the so-called Trotter formula 



+ oir') , (4) 



where r is the time step. For a time interval t — rn, one can combine more 
efficiently n — 1 half steps and obtains 

^t{A+B) ^ ^rA/2 ^^rB^rA^ (""D ^rB^rA/2 (g) 

which essentially corresponds to the leap-frog (or Verlet) algorithm. 
Let us now consider a generic reaction diffusion equation 

V = G(V) + DV^V (6) 

where V = V(x, t) is a classical field. The most used techniques to inte- 
grate such a partial diferential equation (PDE) are the so-called time-splitting 
pseudo-spectral codes llil. These schemes are among the most stable and effi- 
cient and we will illustrate it in the case of the leap-frog time splitting 
Referring to (||) , the operators are identified in the following manner 

AY = G(V) , BY = DV^Y (7) 

where now A\s a non-Hnear local (in space) operator and B a linear non-local 
one. One should now solve separately the evolution equation involving the 
operator A and B and then apply formula (||) . The solution of the non-linear 
part is usually performed by standard integration techniques for ODE, but in 
some cases it can be even solved exactly (as it will be shown for the CGLE). 
To integrate the part involving the spatial derivatives, we consider the Fourier 
transform of the field V(p, t) and of the corresponding evolution equation 

V = -p^DY (8) 

whose solution is simply 

V(p, t + r)^ exp [V^t] V(p, t) (9) 
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In order to obtain the solution in the real space it is now sufficient to back Fourier 
transform the signal. Therefore, an integration over a time At with discrete time 
step T = Ai/n and adopting a leap-frog scheme can be summarized as 



V(r,t + At) = <^e* 



p-l^-p'Drp^rA 



(n-1) 



p~l^-p Drp^iA^Y{^^t) (10) 



where F indicates the Fourier transform and the reverse transform. 
The integration method we propose follows exactly the same scheme as the one 
just shown, but the step involving spatial derivatives is directly solved in the real 
space. As previously noticed, integration in Fourier space reduces to a simple 
product (see (^); in the direct space this corresponds to perform a convolution 
integral of the signal with an appropriate kernel 

V{r,t + T)= J dsK{s,T)V{r~s,t) , (11) 

where r = (xi, x^) in the three dimensional case. The expression of the 
kernel in 3d for the reaction-diffusion systems is straightforward 



ADt 



(12) 



In order to implement this scheme on a computer, we have to deal with a 
discrete time step t and with a discrete spatial resolution A = (Ai,A2,A3). 
The discrete version of the convolution integral (O) reads 



V(l,z + 1)= J2 ^(n)V(l-n,z) (13) 

where r = (Aini, A2n2, A3713) and n — (ni, 712,713) is the spatial index, while i 
is the temporal one {t = ir) and the sum runs over all the lattice. This approach 
maintains exactly the same accuracy of the spectral codes, but it is inefficient in 
terms of CPU time. Since the kernel decays rapidly along any spatial direction. 



a first naive improvement could consist in evaluating the sum in ( |13| ) only over 
q nearest neighbours of each considered site, but this approach, for precisions 
comparable to those of pseudo-spectral codes, is also quite time consuming 0. 
A better strategy consists in truncating the sum in Eq. ( |l3|) and in substituting 
the discretized kernel with a set of unknown coefficients C(n), that should be 
conveniently determined Ii3 

V(l,z + 1)= C(n)V(l-n,z) (14) 

m ,712 ,n3 = — Na 

where Nc is the number of convolution channels considered along each spatial 
direction. We have devised the following method to determine the {2Nc + 1)'^ 
coefficients C(n) for a d-dimensional lattice: 
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• let US first make the ansatz that the field can be well approximated, on the 
considered grid, by decomposing it on {2Nc + l)'^ elements of the following 
Fourier basis 



exp 



J 

L h 



hO-hXh 



where kh = —Nc, . . . ,Nc and au is a set of parameters to be determined. 
On this basis the field can be expressed as 



V(l,z)= ^ V({fc,},i) exp 



h=l 



• the time evolution of each Fourier components is determined by the time 
propagator dsh 



exp ^-'^{khahyOi 



• substituting in eq. ( p. 31) the field decomposed on the Fourier basis and 
imposing that the kernel should give the exact time evolution on these 
modes, one is left with the following set of equations : 



E 



exp 



-j E khahrihAh 



Mkh 



h=l 



C{n) = exp [-(Y^ikhahfDT^ (15) 



from which one can obtain the coefficients C(n). 



In isotropic problems the spatial resolution does not depend on the consid- 
ered axis (i.e. — A V/i) and in that case symmetry reasons suggest that 
ah — aih. Moreover, the symmetry of the kernel allows to reduce the {2Nc + \-)''' 
system to a system of {{Nc + d)\/[Nc\d}] — 1} equations. This because the ele- 
ments of the kernel are related by the following symmetry relations 



C(ni,rt2,n3) = C'{[n\ 
and by the normalizations condition 



nl]) 



(16) 



C(0,0,0) = 1- C'(ni,n2,«3) (17) 

m, "2, "3^^(0,0,0) 

However, the problem of the choice of the parameters {at} is left unsolved. Its 
values are determined by demanding minimization of the error due to the use 
of a Fourier basis that is not the complete one. The a^-parameters should be 
determined for each choice of the time step, of the spatial resolution and of the 
diffusion coefficient. An empirical strategy to choose the optimal a^-values, that 
turns out to be essentially correct, consists in finding the values that minimize 
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the quadratic sum SSQ = J2ni n2 n3=-N^ |C('^i, ?i2, «3)P, with the additional 
constraint that the first 6 cumulants of the discretized kernel should coincide 
(within an arbitrary precision) with those of the true kernel (|l2|). In Fig. |l| it is 
shown that such empirical recipe indeed coincides with the request of a minimal 
integration error. 

2e-04 I ■ . ■ . ■ 1 




-1e-04 I ' ' ' ' ' 1 

0.1 0.3 0.5 0.7 

a 

Figure 1: Dependence of the integration error (dashed line) and of the SSQ 
(reported in arbitrary units) (circles) on the choice of the parameter a for the 
one-dimensional FHNE, with A^^ = 0.5, t = 0.10, iVc = 5 and no-flux boundary 
conditions. The differences between the values of the 4*'*- (triangles) and 6*''- 
momenta (squares) of the corresponding kernel and those of the true one (^ 
are also reported. Numerical instabilities have been usually found for small 
a- values, that in the present case correspond to a < 0.08. 

3. Specific applications 

In this Section we will describe in detail the implementation of our integration 
scheme for the CGLE and the FHNE, both in one and two dimensions. 
Let us first consider the FHNE, in this case we have chosen values of the pa- 
rameter for which the system is in an excitable state, namely a — 0.015, = 1, 
e = 0.006 and D = 1. We have performed the integration of the nonhnear part 
with a simple Euler scheme, while for what concerns the diffusive part of the 
equation we have simply to integrate the equation: 

ut = D\/^u (18) 

This can be done by adopting our scheme where the kernel is given by the real 
function reported in (^2|). In Fig. || a comparison between the "real" kernel 
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K{r) and its optimized expression C(k) on the spatio-temporal grid, obtained 
by solving the system ([l^), is reported. 

A more interesting situation is represented by the CGLE, because this time one 
is able to solve exactly the differentiaJLequation involving the nonlinear operator. 
This result was firstly reported in t§, where the authors have shown that, at 
least for the CGLE, the precision of standard time-splitting pseudo-spectral 
integration schemes is superior to that associated with predictor-corrector and 
Runge-Kutta algorithms. In order to solve such ODE one has to write separately 
the following equations for the phase 

'-^^c,Ar,t) (19) 

and for the amplitude 

^''^'^'^ =2[p\r,t)-p*ir,t)] . (20) 



dt 

The solution of Eq. (|20|) is simply given by 



p(r, r) = l/v/c-2-(l/p2(r,0)-l) + l (21) 

and inserting this into eq. one obtains 

^(r,r) =7A(r,0)+C3{r + ln[p(r,0)Mr,r)]} . (22) 

For what concerns the integration of the diffusive part, it should be noticed 
that this time the kernel is complex and it has the expression 

Kir) - (iy^\-''^M' [cos(7^|rn - jsin(7,|r|2)] (23) 

where 7 = 7r - hi 1/[4t(1 + jci)]. 



4. Treatment of boundary conditions 

The most commonly used boundary conditions are: periodic (P), no-flux 
(NF) (where the first derivative of the field vanishes at the boundaries) and fixed 
(FO) (where the field should have a fixed value at the boundaries, typically set 
to zero). To treat these different cases with pseudo-spectral codes, one is obliged 
to use different Fourier basis llj : namely, the complex basis for P, the cosine 
basis for NF and the sine basis for FO. Instead, with the present algorithm 
the same kernel can be used for each of the above cited cases and different 
boundary conditions are treated just manipulating in different ways the field 
values at boundaries. This allows to treat, e.g. a discretized two-dimensional 
field with different boundary conditions on each border or even with various 
boundary conditions on each single border. 
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Let US consider a classical field {V(fc, i)}, with A; = 0, . . . , iV — 1, discretized 
on a one-dimensional grid of resolution Ax at a fixed time i. Problems arise in 
the evaluation of the convolution ( p^ when k > N — Nc — 1 and k < Nc, when 
points "outside" the chain should be considered. For P conditions the field is 
replicated out of the chain following the simple rules 'V{N + 'm, i) = V(m, i) and 
V(— m,?) = V(iV — m,i). In such case it should be noticed that the first point 
of the grid is located at any point of the chain, while the last point is placed at 
a distance L, where N x — L. The NF conditions are treated differently, in 
this case we have Y{N + m,i) = V{N — m + 1, i) and V(— to, i) = V(to + 1, z), 
moreover the first and the last point of the grid are located inside the chain at a 
distance {Ax)/2 from its extrema. Finally, the PO conditions are implemented 
as follows : V{N + to,?) = — V(iV — m,i) for m ^ and V(7V, i) — 0, while 
V(— TO,i) = — V(m, i). In this last situation, the first point of the grid is 
located at the beginning of the chain, while the last point is within the chain at 
a distance A^^ from the extremum. 

5. Precision and Performances of the Code 

In the present Section we will compare our algorithm with standard time 
splitting codes Efl for the one- and two-dimensional FHNE and CGLE. It should 
be stressed that-the integration algorithm introduced here and the usual time 
splitting codes Il3 differ only in the treatment of the spatial derivatives. There- 
fore for what concerns the analysis of the precision and of the performances we 
will concentrate on such integration step. 

In order to give a measure of the precision of the considered algorithms, we 
have estimated the mean square deviation (MSD) between an orbit obtained by 
the algorithm under test and a reference orbit, that is assumed to be "exact". 
In particular, the MSD has been evaluated over a fixed time interval t and 
averaged over a total integration time T — nr. Once the test orbit has been 
periodically synchronized at the reference one at times t = tot (to = 1, . . . , n). 

In the considered time-split integration schemes we can identify two different 
type of errors, one due to the spatial resolution and the other to the temporal 
one. The second kind of error can be measured considering a test and reference 
orbit obtained with the same spatial resolution but with different time integra- 
tion steps T (obviously that of the reference orbit should be taken much smaller, 
typically we have considered ~ t/5). The "temporal" MSD are obviously identi- 
cal for our algorithm and for the standard pseudo-spectral ones. Moreover, with 
the spatial resolution considered (that corresponds to typical values employed 
in literature) this error is dominant with respect to the "spatial" one. Data are 
reported in Table I. 

For what concerns the MSD associated to the spatial resolution, it is clear 
from the data reported in Table II, III, IV and V that this error is different 
for the 2 algorithms and it depends strongly on the number of convolution 
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T 


Id 


2d 


0.1 


.164E-03 


.742E-02 


0.05 


.741E-04 


.301E-03 



Table 1: MSD associated to the temporal integration for the FHNE: the data 
refers to a spatial resolution = 0.5 and no-flux boundary conditions. 

cliannols TV,; (whon our algorithm is considered). Our algorithm being dovised 
in order to reproduce the precision of pseudo-spectral codes (that we will term 
FFT), one expects that the MSD associated to our algorithm will be bigger 
than that associated to the FFT code for any value of Nc. However, already for 
Nc — 20 — 30 in the one- and two-dimensional cases the precision of the FFT 
code is recovered. 



A^ 


1.0 


2.0 


iVe = 5 


2.961E-07 


9.399E-05 


A^c = 10 


7.770E-08 


8.570E-05 


Nc = 20 


6.593E-08 


7.289E-05 


FFT 


6.727E-08 


7.513E-05 



Table 2: MSD associated to the spatial integration relative to our algorithm for 
various values of Nc, the last row refers to usual pseudo-spectral code. The data 
have been obtained for the one-dimensional FHNE with time step r = 0.05 and 
no-flux boundary conditions. 



Ax 




1.0 


2.0 


Nc = 


3 


9.115E-07 


7.619E-05 


Nc = 


5 


1.256E-07 


4.913E-05 


Nc = 


6 


6.884E-08 


4.264E-05 


Nc = 


7 


5.128E-08 


4.292E-05 


Nc = 


8 


3.596E-08 


4.066E-05 


Nc = 


10 


2.838E-08 


4.102E-05 


FFT 


1.942E-08 


2.835E-05 



Table 3: MSD associated to the spatial integration relative to our algorithm for 
various values of Nc, the last row refers to usual pseudo-spectral code. The data 
have been obtained for the two-dimensional FHNE with time step r = 0.10 and 
no-flux boundary conditions. 

The CPU time required by the pseudo-spectral code will increase with the 
linear dimension N (where the total number of mesh points in d dimension is 
N'^) as N'^ln{N), instead our algorithm will scale as N'''f(Nc). In particular 
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2.00 


1.00 


0.50 


Nc = 


: 6 


0.125E-04 


0.266E-06 


0.123E-08 




10 


0.955E-05 


0.187E-07 


0.120E-11 




20 


0.227E-05 


0.109E-08 




Nc^ 


30 


0.211E-05 


0.617E-09 


0.109E-13 


Nc^ 


40 


0.455E-05 


0.525E-09 


0.921E-14 


Nc = 


50 


0.212E-05 


0.521E-09 


0.866E-14 


FFT 


0.212E-05 


0.521E-09 


0.866E-14 



Table 4: MSD associated to the spatial integration relative to our algorithm 
for various values of Nc, the last row refers to usual pseudo-spectral code. The 
data have been obtained for the one-dimensional CGLE with time step r — 0.05, 
periodic bounday conditions and for parameter values ci = 3.5 and C3 = 0.9. In 
this case the temporal-MSD is 0.161i? — 02. The data of this Table have been 
previously reported in llJ. 



0.25 0.50 0.75 



Nc = 5 
Nc = 10 
Nc = 15 
Nc = 20 
FFT 



0.654E-06 
0.265E-07 
0.113E-10 
0.254E-13 



0.510E-04 
0.463E-05 
0.359E-06 
0.308E-06 
0.225E-06 
0.214E-06 



0.316E-03 
0.724E-04 
0.264E-04 
0.239E-04 
0.724E-04 
0.256E-04 



Table 5: MSD associated to the spatial integration relative to our algorithm 
for various values of Nc, the last row refers to usual pseudo-spectral code. The 
data have been obtained for the two-dimensional CGLE with time step r = 0.05, 
periodic boundary conditions and for parameter values ci — 3.5 and C3 — 0.9. 
In this case the temporal MSD is 0.502£; - 03. 

for the FHNE in the one-dimensional case the ratio between the CPU time 
requested by the code introduced here an that associated to a standard FFT 
routine is 77Vc/ln(iV), where 7 ~ 0.15. Therefore for typical values {Nc = 10, 
N = 2048) such ratio is 0.2. In the two-dimensional case, the ratio can be 
expressed as 7[(7Vc + l)(A^c + 2)/2- 1]/ ln(7V) and for the FHNE we have 7 - 0.3. 
In this last situation for A^c = 4 and N — 2048, the CPU time ratio has a value 
0.6. Thus our algorithm is noticeably faster than usual pseudo-spectral codes 
in one dimension, while in two dimension it becomes convenient to use for grids 
larger than N - 2048. 

6. Conclusions 

In the present paper, an innovative integration scheme for reaction-diffusion 
PDEs has been reported. We have shown that such scheme is competitive, both 
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in one and two dimensions, with respect to usual time-splitting pseudo-spectral 
codes. It combines high levels of accuracy with reduced cost in terms of CPU 
time and an extreme portability, being applicable with a quite limited effort 
to systems with different boundary conditions. Future improvements of the 
present algorithm should consist in determining an automatic procedure for the 
choice of the a^-parameters. We plan to implement this algorithm in a parallel 
enviroment and to extend its applicability also to more general classes of PDEs. 

Acknowledgements 

The first version of this integration scheme has been developed by one of us, 
for the one-dimensional CGLE, mainly JbUDwing Peter Grassberger's ideas and 
in collaboration with Helge Frauenkron ll3'tZI. We wish also to thank Markus Bar 
and Roberto Genesio for useful discussions. A.T. acknowledges the hospitality 
of family Frese in Wuppertal (Germany) during the final write up of this text 
and Leonardo Torcini for providing him a lap-top. 

References 

1. E. Bettel heim and B. Lehma nn, Ann. Rev. of Comp. Phys., 7 to appear in 1999 
(preprint |3ond-mat/9904205| ). 

2. M.C. Cross and P.C. Hohenberg, Rev. Mod. Phys. 65 851 (1993). 

3. P. Grinrod, Patterns and Waves. The theory and application of reaction- 
diffusion Equations, (Clarendon Press, Oxford, 1991). 

4. W. Van Sarloos, Phys. Rep. 301 9 (1998). 

5. J. Tyson and J.P. Keener, Physica D 32 327 (1988). 

6. A.T. Winfree. Chaos 1 303 (1991). 

7. Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, (Springer Verlag, 
Berlin, 1984). 

8. H. Sakaguchi, Prog. Theor. Phys. 84 792 (1990); B.J. Shraiman et al, Physica D 
57 241 (1992); H. Chate, Nonlinearity 7 185 (1994); D.A. Egolf and H.S. Greenside, 
Phys. Rev. Lett. 74 1751 (1995); P. Manneville and H. Chate, Physica D 96 30 
(1996); R. Montagne et al, Phys. Rev. Lett. 77 267 (1997); A. Torcini, Phys. Rev. 
Lett. 77 1047 (1997). 

9. J.D. Murray, Mathematical Biology, (Springer Verlag, Berlin 1993). 

10. L. Verlet, Phys. Rev. 159 98 (1967). 

11. M.P. Allen and D.J. Tildesley, Computer simulations liquids, (Clarendon Press, 
Oxford, 1987). 

12. H. Yoshida, Phys. Lett. A 150 262 (1990); E. Forest and R.D. Ruth, Physica D 
43 105 (1990); J. Candy and W. Rozmus, J. Comp. Phys. 92 230 (1991); J.M. 
Sanz-Serna, Physica D 60 293 (1992); M. Suzuki, Int. J. Mod. Phys. C 7 355 
(1996). 

13. R.I. McLAchlan and P. Atela, Nonlinearity 5 541 (1992). 

14. W.H. Press et al.. Numerical Recipes, (Cambridge University Press, Cambridge, 
1986). 

15. H. Frauenkron and P. Grassberger, Int. J. Mod. Phys. C5 37 (1994). 

16. A. Torcini, H. Frauenkron, and P. Grassberger, Physica D 103 605 (1997). 

17. A. Torcini, H. Frauenkron, and P. Grassberger, Phys. Rev. E 55 5073 (1997). 



M. mm, A. Torcini & S. Ruffo 
D. Goldman and L. Sirovich, Quart. Appl. Math. 53 315 (1995). 



An Integration Scheme for Reaction- Diffusion Models 13 





Figure 2: Comparison of the true kernel with the optimized one for three differ- 
ent values of Nc = 3, 5, 10 for the one-dimensional FHNE, obtained employing 
the following parameters r = 0.05, A^; = 0.25 and a = 0.10. 



